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ABSTRACT 


A  new  method  is  presented  for  fitting  polynomial 
splines  to  n  equispaced  data.  Using  Jain's  [32]  cyclic 
decomposition  of  banded  Toeplitz  matrices,  we  show  that  the 
operations  can  be  performed  by  n-point  Fast  Fourier 
Transforms  ( FFT ) .  Thus,  the  use  of  parallel  processing  FFT 
techniques  provides  a  speed  of  O^log^n),  independently  of 
the  degree  of  the  spline.  Explicit  solutions  are  derived 
for  the  cubic,  quartic  and  quintic  spline. 


Research  supported  by  the  Office  of  Naval  Research  through 
Grant  N00014-82-K-0409 . 


Introduction 


In  many  applications  there  appears  the  need  to  represent 
a  set  of  raw  data  by  fitting  to  them  a  smooth  function  or  set  of 
functions.  A  great  amount  of  theoretical  work  has  been  done  within 
the  framework  of  approximation  theory  in  studying  properties  of 
curve  fitting  for  various  families  of  approximating  functions. 

One  of  the  most  attractive  and  well  structured  families  of 


approximating  functions  are  the  splines.  They  have  been 
extensively  studied  in  the  mathematical  literature,  for  example 
in  [1]  -  [7]  and  elsewhere. 

Splines  have  been  very  useful  in  statistics.  References 
[7]  -  [18]  represent  the  most  significant  spline  application 
papers  in  the  statistical  literature. 

" —  In  the  engineering  literature,  spline  functions  have  been 
used  as  approximating  tools,  in  the  areas  of  Systems  (.[19]  -  [22]) 
and  Pattern  Recognition.  ([23]  -  [30]). 

In  the  present-  paper  we  concentrate  on  a  fast,  parallel 
computation  technique  for  fitting  a  spline  to  n  equispaced  data 
points.  Existing  techniques  can  fit  splines  in  0 (n)  time,  with 
recursive  (serial)  processing  of  the  data,  as  in  [16].  The 
new  technique  is  based  on  the  use  of  Fast  Fourier  Transform,  and 


its  parallel  processing  capability.  As  a  result,  our  technique 

t  a. 

achieves  the  fitting  of  a  spline  in  Odogjn)  time.  The  n 
data  points  must  be  equispaced. 


Let  ((x^,y^)  i  =  0,  1,  . ..,  n  }  be  a  set  of  data  points 


with  a  =  xrt  <  x,  <  ...  <  x  =  b.  We  would  like  to  fit  to  the 

u  i  n 


data  a  function  S(x)  that  has  two  derivatives.  We  require  that 
S(x^)  =  y^,  i  =  0,  1,  ...»  n  and  that  under  the  above  requirement 
S  minimizes  the  integral 


I2(S)  =  /  [S" (x) ]  ax- 


(1) 


under  the  conditions  S' (a)  =  S'{b)  =  0. 

In  the  theory  of  spline  functions  it  is  shown  [2]  that  the  above 
constrained  minimum  is  achieved  when  S (x)  is  a  set  of  piecewise 
cubic  polynomials  with  continuous  first  and  second  derivatives 


at  the  points  {x^  i=l,  ...»  n=l}.  We  will  be  concerned  here  only 


with  uniformly  spaced  x^'s,  hence  x^  =  xq  +  iht  i  =  0f  1,  . . . ,  n. 


h  =  increment.  The  continuity  requirement  demands  that; 


S' (x. )  =  S' (xt) ,  S"(xi)  =  S"(xt) 


(2) 


for  i  =  1,  2,  . . . ,  n-1, 


We  denote  M.  the  second  derivative  of  S(x)  at  x. .  Then 

i  l 


S"(x)  =  Mi_i  (x^-x) h-1  +  (x-xi_1)h_1,  xi_i  1  x  xi 


(3) 


Integrating  (3)  twice  and  using  the  conditions  (S(Xj.)  =  yi }  we  obtain; 


S  (x)  =  Mi_1(xi-x)  3(6h)_1  +  Mj.  (x-xi_1)  3(6h)_1  + 


+  (h  1yi  -  hMi6~1)  <x-xi_1)  +  (h  1yi_1  -  hMi_16_1) (x^-x) 


xi_l  i  x  <  xi 


(4) 


-J 
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Differentiating  (4)  once  and  using  the  continuity  of  the  first 
derivative  for  i=l,  n-1,  we  find: 

Mi-1  +  4Mi  +  Mi+l  =  6h~2(yi-i  "  2yi  +  yi+l)  i=1'  2'  •••'  n_1 

(5) 

We  assume  M  and  M  are  known.  Then  the  set  of  unknowns  is 
o  n 

MjM2  ...  Mn-1  ,  and  we  have  a  set  of  n-1  equations  from  (5)  with 
an  equal  number  of  unknowns .  In  matrix  notation  the  set  of 
equations  is : 

TM  =  6h'2H  (6) 

where 
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y0  -  2?!  +  y2  -  h^e'1 

yx  -  2y2  +  y3 

y2  -  2y3  +  y4 


y„-3  "  2yn-2  +  yn-l 

y„-2  -  ^n-l  +  yn  -  A.6'1 


h 


n-1 


and  T  is  an  (n-l)x(n-l)  matrix,  H  is  an  (n-1)  vector.  Clearly  the 
basic  problem  here  is  the  efficient  inversion  of  T  and  the 
multiplication  T-1H.  We  define  the  following  (n-1) x (n-1)  matrix: 


1  P 

p  1  p  0 
pip 

T  =  *  ‘ 

p  .  . 


1  P 
P  1 


It  is  a  known  result  from  [31]  that  the  eigenvalues  qk  and  eigen¬ 
vectors  V.  of  T  are : 

K  p 

qk  =  1  +  2p  cos  (kirn-1)  (8) 

=  [sin(kwn  ^)  ,  sin(2k7rn  ^),  sin  ( (n-1)  kirn  ^)]/2/n 

k  =  1 ,  . . . ,  n-1  (9) 

If  we  divide  T  by  4,  we  have  the  special  case  of  (7)  with  p  =  1/4. 
Using  the  eigenvector-eigenvalue  expansion  of  T,  we  have: 
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T  =  4  l  [1  4-  0.5  cos  (kirn-1)  ]  V.  v.  (10) 

k=l 

The  inverse  of  T  has  the  same  eigenvectors  and  inverse  eigenvalues. 
Therefore, 

T*1  =  1/4  Y  [1  +  0.5  cos(kTTn_1)  rV^V.  (11) 

k=l  K  K 


The  equation  for  the  coefficient  vector  M  is: 


M  =  1.5h-2  s.  V. 


n-1 


k=l 


k  k  k 


(12) 


where 


=  [1  +  0.5  cos  (kirn  ^)  ]  "*■ 

The  matrix  has  elements: 

{2n  1sin(kmTrn  *)  ,  sinfkqTTn  \  m,  q=l,  ...,  n-1) 


(13) 


(14) 


If  we  define  F  =  (f.  ...  f  .), 

l  n-i 

n-1 

1 

q=l 

then  the  mth  component  of  the  vector  M  is 


=  n-1/2 

Lk  _i_,  q 


f,_  =  n  “  J  h  sin(kqirn  1)  ,  k=l,  . ..,  n-1 


(15) 


M  =  3h  2n  ^2  y  s. f.  sin(km7rn  ^)  (16) 

m  k=l  K  K 

Computations  (15)  and  (16)  are  slightly  modified  versions  of 
Finite  Fourier  Transforms,  and  each  of  them  can  be  completed  in 
0(log2n)  time,  if  Fast  Fourier  Transform  methods  are  used.  In 
other  words,  F  is  produced  from  H  by  a  Fast  Fourier  Transform,  and 
then  M  is  produced  by  the  vector  (s^f^,  .  ..,  sn-i^n-i^  by  another 
Fast  Fourier  Transform.  Hence  the  computation  of  M  requires 
computation  time  of  the  order  0(log2r»). 
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The  reason  for  obtaining  such  a  simple  solution  is  the  fortunate 
event  of  the  eigenvalues  and  eigenfunctions  of  T  being  expressed 
in  a  closed  form  that  relates  them  to  the  Fast  Fourier  Transform. 

For  higher  order  splines  our  solution  requires  to  develop 
methods  based  on  the  approximation  of  Toeplitz  matrices  by  circulant 
ones,  to  be  presented  in  the  next  section.  We  will  utilize  a 
method  of  partitioning  and  cyclical  decomposition  of  banded  Toeplitz 
matrices,  which  is  due  to  A.K.  Jain  [32]. 


Circulant  and  Toeplitz  Matrices 
Part  of  the  exposition  in  the  present  section  follows  Gray 

[33]  . 

A  circulant  matrix  C.  is  one  having  the  form: 


C0  C1  c2  ...  .  cn-1 

cn-l  c0  C1  c2’  ’  *  cn-2 

cn-l 


(17) 


'n-1 


c2 

C1 

C0 


The  eigenvalues  qm  and  eigenvectors  of  C  are  the  solutions  of 


CV  =  qV,  V  =  [vQ,  ...,  vn_1]t 


(18) 


or  equivalently  of  the  n  difference  equations 
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m-1  n-1 

Jn  cn-m+k  vk  +  l  °k-mvk  “  <Jvm'  m=0-  1 . "-1  1191 

k=0  k=m 

It  is  easily  verified  for  any  m=0 ,  1,  ...»  n-1,  that  =  exp{-2irimkn  "*■} 
is  a  solution  to  (19) ,  resulting  in  the  eigenvalues 


q  =  7  c  exp{-2iriinkn  } 

m  k==0  * 

(i  =  /=!) 


and  the  corresponding  eigenvectors 


V  =  n  ^2[1,  exp  (-2-irimn  ^)  , 
m 


,  exp  (-2Trim(n-l)  n  *)  ] (21) 


We  can  now  write 

-l  n  ~  ^  *+- 

C  =  n  I  q  V  V  ' 
m  m  m 

m=0 


c-1  =  n"1  y  q“Lvv* fc 
a  in  m  m 
m=u 


We  observe  that  all  circulant  matrices  have  the  same  set  of 
eigenvectors.  Also,  the  inverse  of  a  circulant  is  also  a 
circulant.  The  multiplication  of  C  1  to  any  vector  H  =  (hg...hn_1)t 
can  be  done  by  Fast  Fourier  Transform  techniques  as  follows: 

Let  Z  —  (zcz1  ...  , 

Z  =  C_1H  (24) 

—  1/9  ^  —  1 

Let  f,  =  n  '  y  h  exp(-2wiskn  x)  (25) 

K  s=0  s 

k  =  0 ,  . . . ,  n-1 

be  the  Finite  Fourier  Transform  coefficients  of  H.  Then,  the 


components  of  Z, 


(26) 


z 


m 


n-1/2 


-1  -1 
I  «3v.  f.exp  ( 2  wimkn  ) 

k=l  K  K 


m  -  0,  ...»  n-1 

are  Finite  Fourier  Transform  coefficients.  Hence,  Z  is  computable 
from  (26)  in  0(log2n)  computation  time.  A  Toeplitz  matrix  of 
order  (n,  m,  s) ,  with  s  <  n,  m  <  n  is  defined  as  an  nxn  matrix  with 
entries  t(k,  j)  such  that 


t (k,  j) 


t(k-j)  for  -s  <_  k-j  m 
0  otherwise 


and  t(m)  ^  0,  t(-s)  5*  0. 


t(0)  t(-l)  .  .  t(-s) 

0 

t(0) 


1 27) 


T 

n 


t  (m) 


t(0)  .  .  t(-s) 


t(m)  .  .  t(0) 


(28) 


With  the  exception  of  the  upper  right  and  lower  left  corners 
Tn  looks  like  a  criculant  matrix;  i.e.,  each  row  is  the  row  above 
it  shifted  to  the  right  one  place.  If  we  fill  in  the  upper  right 
and  lower  left  corners  by  the  appropriate  entries,  we  can  make  Tn 
exactly  a  circulant.  Define  the  circulant  matrix  C  in  this  way 


c 


k 


t (-k)  ,  k  =  0,  . 
t(n-k) ,  k  =  n-m, 


,  s 
n-1 


(29) 


0 


otherwise 
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C,  defined  as  above,  is  a  prime  candidate  for  approximating 
T  .  The  rationale  for  such  an  approximation  is  the  desire  to 
approximate  any  operation  Tn^Y,  Y  =  vector,  by  the  operation  C^Y, 
which  can  be  performed  in  0(lognn)  computation  time  using  the  Fast  Fourier 
Transform. 

Let  D  be  the  difference 

D  =  C-T  (30) 

n 

D  has  the  form 


D 


0 

P 


Q 

0 


where  Q 


t(m)  t(m-l) 
t  (m) 

0 


t(l) 
t (m-1) 

t  (m) 


t(-s)  0 

t(l-s)  . 

t(-l)  .  .  t(l-s)  t(-s) 


Now  suppose  that  we  want  to  solve  the  equation  TnX  =  Y,  where 

X,  Y  are  n-vectors.  Substitute  T  =  C-D.  Then  we  have 

n 

CX  =  DX  +  Y  or  X  =  C_1DX  +  C-1Y  (31) 

We  partition  X  and  C  ^Y  as  follows: 


*1 

~Wl" 

-1 

X2 

,  C  XY  = 

W2 

X3 

W3 

—  — 

_  _ 

(10) 


The  method  of  solution  of  equations  (34)  -  (36)  through 
partitioning  and  cyclical  decomposition  of  the  banded  Toeplitz 
matrices,  is  due  to  A.K.  Jain  [32]  and  is  a  very  efficient  method 
for  solving  a  system  of  linear  equations  when  the  matrix  Tn  is  of 
the  banded  Toeplitz  form. 

Using  parallel  processing  techniques  through  Fast  Fourier 
Transform  architecture,  the  solution  of  equations  (34)  -  (37)  is 
achieved  in  (s+m) ^  +  0(21og2n)  time. 


i 

( 


(11) 


In  the  present  section  we  wil.  use  the  method  of  approximating 

a  Toeplitz  matrix  by  a  circulant  to  solve  efficiently  the  spline 

fit  problem  for  higher  degree  splines.  Explicit  solutions  will 

be  given  the  quartic  and  quintic  spline,  which  are  the  simplest 

and  hence  more  useful  higher  degree  splines.  To  avoid 

proliferation  of  notation,  we  will  use  to  denote  the  kth 

derivative  at  the  knot  x^  of  a  k+1  degree  spline.  Hence,  k=2,  3,  4 

correspond  to  the  cubic,  quartic  and  quintic  spline.  The  number 

of  knots  used  will  be  n+k  for  the  k+1  degree  spline,  so  that  k 

of  them  serve  as  boundary  points  with  parameters  assumed  known, 

leaving  exactly  n  unknown  parameters.  They  form  the  vector 

M  =  [M.  ...  M  ]fc.  Hence  the  matrix  to  be  inverted  will  be 
I  n 

always  nxn.  The  vector  M  of  kth  derivatives  together  with  the 
continuity  of  the  first  k  derivatives  are  sufficient  to  define 
the  spline  completely. 

We  consider  now  the  quartic  spline.  Let 


I4  =  {Xj  =  jh,  a  =  -h,  b  =  (n+l)h,  j  =  -1,  0,  1,  ...,  n,  n+1} 

(38) 

be  a  set  of  n+3  equispaced  knots  on  [a,  b] .  A  quartic  spline  S(x) 
defined  on  [a,  b]  is  a  piecewise  fourth  degree  polynomial  between 
knots  (x.x...)  that  fits  a  set  of  data  (S(x.)  =  y.,  jel.}  and  has 
continuous  the  first  three  derivatives  at  the  knots,  i.e. 

S(k)(x7)  =  S(k)(Xj),  k  =  0,  1,  2,  3,  jcl4  (39) 

Let  M.  =  S (3) (x~)  =  S (3) (Xj) ,  jel4 


(40) 


(12) 


We  have  (n+3)  parameters  {  M ^ ,  jel^  }.  The  quartic  spline 
requires  specification  of  the  following  3  boundary  numbers: 

S' (a) ,  S" (a) ,  S' (b) , 

It  is  shown  in  [4]  that  the  3  boundary  conditions  specify  uniquely 
the  "boundary  parameters"  M^,  MQ,  Mfi+1.  Therefore,  we  will 
consider  them  known.  It  is  further  shown  in  [4]  that  M ^  satisfy 
the  following  equations : 

”k-2  +  11Mk-i  +  11Mk  +  Mk+i  a  24h'3(_yk-2  +  3yk-i  ■  3yk  +  yk+i} 
for  k  =  1,  2,  . . . ,  n  (41) 

We  have  now  n  equations  and  an  equal  number  of  unknown 
parameters,  constituting  the  vector  M  =  [Mx  ...  M  ]fc. 

We  define  the  n  constants 

-h3MQ.  11/24  -  h3M_x/24  +  (y2  -  3yx  +  3yQ  -  y^)  ,  k=l 
-h3MQ/24  +  (y 3  -  3y2  +  3yL  -  yQ) ,  k=2 

: 

(yj+1  -  3Yj  +  3yj_1  -  y^_2),  3  <  k  =  j  <  n-1 

h\+l/24  +  (yn+l  "  3yn  +  3yn-l  "  yn-2>  '  k=n 

Let  F  =  (f.  .  .  .  f  ) t 
i  n 

We  also  define  the  nxn  Toeplitz  matrix 


(42) 


(13) 


(43) 


Then,  the  vector  M  from  the  solution  of  equations  (43)  is: 

M  =  24h-3T_1F  (44) 

At  this  point,  we  are  ready  to  apply  directly  the  procedure  of 

Section  II  because  T  is  banded  Toeplitz.  If  we  identify  T  with 

T  from  eq.  (30),  we  have:  s  =  1,  m  =  2,  t(0)  =  t(l)  =  b, 
n 

t (2)  =  1,  t(-l)  =  1.  The  matrices  P,  Q  are: 

Q  *  J  ^  ,  P  *  1  (scalar)  (45) 

Also, 

Y  =  24h  3F,  x^  =  scalar 

x^  =  two  dimensional  column  vector.  Here  the  unknown  vector  M  is 
identified  with  X.  Using  equations  (34)  -  (37),  the  solution  for 

.  3 

M=X  is  immediate.  The  time  required  for  the  operations  is  3  +0(2  log2n) 
Finally,  we  will  derive  the  solution  to  the  quintic  spline 
fit  problem.  Let 


I,.  =  {Xj  =  jh,  j  =  -1,  0,  1,  ...,  n,  n+1,  n+2,  a  =  -h,  b  =  (n+2)h) 

(46) 

be  a  set  of  (n+4)  equispaced  knots  in  [a,  b] .  A  quintic  spline 

S(x)  defined  on  [a,  b]  is  a  function  that  has  continuous  the  three  first 

derivatives  on  la,  b],  fits  a  set  of  data  y^  =  Sly^),  jel5,  and 


(14) 


minimizes  the  integral 

As<3>  (x)  ]2dx 
a 

under  4  specified  boundary  conditions: 


(47) 


S' (a) ,  S"  (a) ,  S’  (b) ,  S" (b) . 
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(The  quintic  spline  corresponds  to  a  differential  operator  L=D  ) . 

In  spline  theory  [4]  it  is  shown  that  the  solution  to  the  above 
constrained  minimization  problem  is  a  piecewise  fifth  order  poly¬ 
nomial  with  the  first  four  derivatives  continuous  at  the  joints. 
Hence, 

SUj)  =  y  j ,  S(k)(x7)  =  S(k)(xt),  k  =  0,  1,  2,  3,  4,  jeI5  (48) 

Let  Mj  =  S(4)(xT)  =  S(4)(x*),  jel5  (49) 

The  4  boundary  conditions  specify  uniquely  the  4  "boundary 
parameters"  M_1,  MQ,  Mr+1,  M  ^  as  shown  in  [4].  Therefore, 
we  will  consider  them  known.  It  is  shown  in  [4]  that  the  following 
set  of  n  equations  is  satisfied  by  the  NT's: 

«k-2  *  26Mk-l  *  66Mk  +  26M*+1  +  “*.2  - 
120h-4tyk_2  -  4yk_x  +  6yk  -  4yk+1  +  yk+2l 

for  k=l,2,...,n  (50) 

Now  we  have  n  equations  and  n  unknown  parameters,  the  components 
of  the  vector  M  =  [M^M2  ...  Mn] .  Let  F  =  [ f ^ f 2  ...  fnl»  where: 


! 


(15) 


(16) 


(54) 


Y  =  120h"4F  (55) 

Using  (39)  we  can  compute  X^,  with  4  3  operations.  Combining 
with  (37)  we  compute  X2-  The  total  computational  time  required  is: 

4 3  +  0 (21og2n) . 

Generally,  to  fit  a  k+1  degree  polynomial  spline  to  n+k 
equispaced  knots  with  data  {y^},  we  have  n+k-1  intervals  and 
n+k-2  interior  points.  For  each  interval,  specification  of  the 
spline  requires  k+2  constants  to  be  solved  for.  Hence  the 
number  of  unknown  parameters  is  (k+2)  (n+k-1) .  The  continuity 
of  (x) ,  j  =  0,  1,  ...,  k  at  the  interior  knots  (joints) 

provides  (k+1) • (n+k-2)  equations.  By  fitting  the  data  points 
{y j }  we  get  (n+k)  equations.  Let  s  be  the  number  of  specified 
boundary  conditions.  In  order  to  have  equal  number  of  unknown 
parameters  and  equations,  we  must  have: 

(k+2) (n+k-1)  =  (k+1) (n+k-2)  +  (n+k)  +  2 

hence 

s=k . 

The  boundary  conditions  can  be  spec  i  fied  either  by  a  number 
of  derivatives  at  the  end  points  a,  b  or  by  a  number  of  boundary 
values  of  the  moments  M j .  For  the  cubic  spline,  k=2,  and  we 
picked  Mq,  Mn+1  as  the  known  boundary  values.  For  the  quartic 


P  = 


1  0 
26  1 


X1  = 


M, 


M, 


Q  = 


1  26 
0  1 


X3  = 


M 


n-1 


n 


(17) 


spline,  we  have  k=3.  The  moments  <M__ ^ ,  MQ,  M  +^)  were  assumed 
known.  For  the  quintic,  we  have  k=4  and  we  used  the  moments 
(M_1#  Mq,  Mn+1,  Mn+2)  as  known  boundary  conditions.  An  alternate 
set  of  boundary  conditions  that  has  been  frequently  used  is  the 
specification  of  a  total  of  k  derivatives  of  S(x)  at  a  and  b. 

There  is  a  one  to  one  dependence  between  the  two  different  boundary 
conditions  [4].  In  the  present  paper,  the  boundary  moments  were 
chosen  as  more  convenient. 

The  kth  derivative  of  a  (k+1)  degree  spline  is  a  piecewise 
linear  function,  expressed  as: 


(x)  =  ML_^h  3  (Xj  -  x)  +  Mjh  ^(x-Xj_^)  for  xj_^  f.  x  i  xj 

(56) 


where 


-  S(k) (Xj) 


Up  to  this  point  we  have  shown  that  the  parameters  (M^J,  which 
are  the  third  derivatives  at  {x^}  for  the  quartic  spline  and  fourth 
derivatives  at  for  the  quintic  spline,  can  be  computed  from 

the  data  y ^  in  parallel  processing  time  of  21og2n  operations. 

Complete  specification  of  the  other  parameters  of  the  spline 
will  be  now  undertaken. 

For  the  quartic  spline,  we  have: 


S(3)(x)  =  M j_^h” 1 (x j  -  x)  +  Mjh-1 (x  -  x..^) 
for  Xj..!  1  x  i  xj 


(57) 


Integrating  three  times,  we  have: 


(18) 


S (x)  =  -(24b)"1  -  X)4  +  (24h)'L  Mj(x  -  xj_1>4  + 

+  2~3hB3  (x-x . )  2  +  h2B?  (x  -  x.)  +  h3B?,  x.  ,  <  x  <  x.  (58) 

3  3  3  3  3  3-1  ~  ~  3 

12  3 

where  (B^,  B^,  B ^ }  are  constants  that  need  to  be  determined  from 

{M^},  (y j } ,  and  the  continuity  requirement  at  the  joints.  We 

take  the  consecutive  derivatives  of  S(x),  write  the  corresponding 

(k) 

expressions  of  S  (x)  for  the  intervals  [Xj_^,  and  [x^ ,  x^+^]  , 

(k) 

then  apply  the  continuity  requirement  for  S  (x)  at  x  =  x ^ , 


for  k  =  1,  2. 

The  following  equations 

result: 

S(2) (x~) 

-  S(2) (xt) ; 

B3l1  =  B3  + 
3  +  1  3 

M. ,  j  =  0,  1, 

...  ,  n  (59) 

S(1)  (xT) 

•  SU,<x+); 

2  2 

B7.  =  B7  + 
3+1  3 

Bj+1'  }  -  0.  1. 

.  . . ,  n  (60) 

The  equations  {S(x^)  =  y^}  provide  the  following  relationships: 

Yj  =  S  (xT ) ;  B 3  =  h-3yj  -  24'V.,  j  -  0,  1,  ...»  n,  n+1  (61) 

S(a+)  -  y_x;  2”3bJ-B2  =  -B3  +  h'3y_1  +  24~1M_1  (62) 

1  2 

We  still  need  to  determine  the  two  initial  values,  (B^,  B^) .  For  this 
purpose,  we  use  one  additional  equation: 


sl*$>  "  v0>  h'3y0  -  -21’1mo  +  2'1bi  "  Bi 

Substituting  (B3,  B2)  from  (59),  (60)  for  j=0, 
we  get: 

2'3bJ  +  B2  =  B3  -  Mq. 13/24  -  h"3y0 

1  2 

From  (62)  and  (64)  we  can  solve  for  (B^,  Bq) . 


in  terms 


Then  B j , 


(63) 


of  B 


O' 


(64) 


are  computable  from  M^. 


(19) 


For  the  quintic  spline,  we  have: 


S(4)  (x)  =  h-1(Xj  -  x)  +  Mjh-1(x  -  x^^) 

for  xj_i  <_  x  <_  Xj 


(65) 


Integrating  four  times,  we  find: 

S(x)  =  (120h)"1  M^tx..  -  x)5  +  (120h)-1  Mj  (x  -  x^)5  + 

-11  1-122  211 
+  6  h  B j (x  -  Xj)J  +  2  h  Bf(x  -  x^ )  +  h  B]7 (x  -  x^) 


+  h4  B4 
3 

for  xj_^  <.  x  <.  xj 


(66) 


12  3  4 

The  constants  now  are  B . ,  B . ,  B ^ ,  B .  .  We  take  the  consecutive 
(k ) 

derivatives  S  (x)  for  the  intervals  [x^_^,  x^l  and  [x^ ,  x^+^] , 
then  apply  the  continuity  requirement  for  S'  '  (x)  at  x  =  x ^ , 
k  =  1,  2,  3.  We  also  apply  the  condition  S(xT)  =  y^.  The  resulting 
equations  are: 


S(x”) 

4 

=  y j ;  b.  = 

h'4yj 

-  120_1M j  j  =  0,  1,  ...,  n+2 

(67) 

S(3>(x-) 

=  s(3)  (xt)  ; 

B'+i 

=  B.  +  M. 

3  3 

(68) 

S(2>(x-) 

=  s<2>(x;); 

Bi*i 

1  2 

=  B  .  .  +  Bv 

3  +  1  3 

(69) 

s(1) (x~) 

=  SU)  (X^)  ; 

Bj+i 

12  -11  -1 
=  B.  +  B7 ..  -  2  L  B .  ,  .  +  12  M. 

3  3+1  3+1  3 

(70) 

for  j  =  0 ,  1 ,  . . . ,  n+1 


12  3 

We  need  now  to  solve  for  the  initial  values  (Bq,  Bq  ,  BQ) , 

which,  together  with  the  equations  (67)  -  (70) ,  will  give  the 

complete  solution.  For  this  purpose,  we  will  formulate  a  set  of 

12  3 

3  equations,  in  which  the  only  unknown  parameters  are  (BQ,  BQ ,  BQ) . 


(20) 


4 

Note  that  B..  are  immediately  found  from  (67).  From  (66)  we  have: 

h"4.S(Xj_1)  =  120-1Mj_1  ~6-1B j  +  2~1B2  -  B2  +  B4  (71) 

for  j  =  0 ,  1 ,  . .  . ,  n+2 

We  apply  (71)  for  j  =  0,  1,  2.  The  result  is: 


y-lh~4 

=  120-1M_1 

~  6_1bJ  +  2_1Bq 

-  B0 

+  B0 

(72) 

yoh"4 

=  120_1M0 

6-1B^  +  2~1B2  - 

B1  + 

B1 

(73) 

ylh~4 

=  120"1M1  - 

6~1B2  +  2~1B2  - 

B2  + 

(74) 

12  3  12  3 

Using  (68)  -  (70),  we  express  (B1#  B^  B1,  B2>  B2 ,  B2)  in  terms 
12  3 

of  (Bq ,  BQ ,  Bq)  and  substitute  them  in  (72)  -  (74).  Thus  we 

12  3 

finally  have  3  equations  with  the  unknown  parameters  (BQ,  BQ ,  BQ) . 
After  the  determination  of  the  initial  values  BQ ,  equations 
(68)  -  (69)  provide  all  the  parameters  B ^ . 

Conclusions 

We  have  presented  an  efficient  algorithm  for  fitting  a 
polynomial  spline  to  a  set  of  n  equidistant  data.  The  algorithm 
exploits  the  parallel  nature  of  Fast  Fourier  Transform  and  the 
form  of  the  difference  equations  relating  the  moments  { M ^ }  to  the 
data  { y j } .  Explicit  formulas  were  obtained  for  the  cubic,  quartic 
and  quintic  spline,  and  the  method  is  readily  extensible  to  higher 
order  and  more  general  splines.  The  method  involves  only  linear 
operations,  and  is  much  simpler  than  previous  ones.  The  time 
complexity  of  the  method  is  0(2  log2n) ,  independently  of  the  order  o 
the  spline. 
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